The Vortex Method With Finite Elements

This work shows that the method of charcteristics is well suited for the numerical solution of first order hyperbolic partial differential equations whose coefficients are approximated by functions piecewise constant on a finite element triangulation of the domain of integration. We apply this method to the numerical solution of Euler's equation and prove convergence when the time step and the mesh size tend to zero. The proof is based upon the results of regularity given by Kato and Wolibner and on L°° estimates for the solution of the Dirichlet problem given by Nitsche. The method obtained belongs to the family of vortex methods usually studied in a finite difference context. Introduction. The vortex method is based on an old concept of fluid mechanics which says that for two-dimensional nonviscous flows the vorticity in the fluid is transported by the flow; thus, if the initial distribution of vorticity consists of a finite number of point vortices, the flow at later times can be found by transport of these point vortices along the streamlines of the flow that they create. In mathematical terms this means that the two-dimensional stream functionvorticity formulation of Euler's equations -^+«Vw = 0, -A* = w, m = VA* in fi X 10, 7T, (1) J 3/ J ' L' u(t = 0) = 2 to,°5(x x,), *|r = *0, where 8 is the Dirac function, T, the boundary of fi, is integrated by (2) ■ This method was first implemented by Christiansen [6] and Chorin [5] and thoroughly tested by Baker [2] on the roll up of vortex sheets. From the theoretical point of view if fi = R2 Hald [10] showed that when (2) is discretized explicitly in time, when the Dirichlet problem is approximated by a suitable discretization of the corresponding Green's function and when the Dirac functions are smoothed by appropriate convolutions then the method converges. In Baker [2] the Dirichlet problem is discretized with finite differences, and as far as we know the convergence is not established in such a case. The present work is based on the rather straightforward observation that the system (2) is perhaps easier to analyze when it is discretized by the finite element Received May 1, 1980; revised June 30, 1980. 1980 Mathematics Subject Classification. Primary 65N30, 65N35; Secondary 76C05. © 1981 American Mathematical Society 0025-571 8/8 l/0000-OO09/$O5.50 119 License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 120 CLAUDE BARDOS, MICHEL BERCOVIER AND OLIVIER PIRONNEAU method than by the two previously mentioned methods because the equation for the characteristics x,(r) can be integrated exactly if VA* is piecewise constant on a triangulation of fi X ]0, T[. However, the error analysis shows that in the finite element context it is no longer feasible to work with Dirac functions; it is better to use a piecewise constant discretization of w°(x). Therefore, we shall not work with point vortices but with a piecewise constant approximation of the vortex field (3) «(*, 0 = 2 «?(,)/(* *,(')) 1 = 1 where I(x x¡(t)) equals 1 if x and x,(r) belong to the same element of the triangulation and zero otherwise, where TV is the total number of elements and where j(i) is the index of the element to which x,(0) belongs. Therefore we will have to compute certain characteristics backward in t in order to define w(x, t) by (3). Thus, although in spirit identical, in practice the present method is substantially different from the point vortex method of [2] and [5]. Both have the advantage of being nondissipative; ours is conservative in a statistical sense only in terms of «°. On the other hand, we do not have to insert new vortices in some regions of the flow as in [6] and the method is more appropriate to smooth flows. But most of all an error analysis will be given and the method is unconditionally stable in time. This, by the way, may also be true of the cloud and cell vortex method [5], [2] as was observed by Baker. The proofs are involved and difficult in their details but the guidelines are simple: we assume that the regularity obtained by Wolibner [14] and Kato [11] for the solution holds. Thus, to measure the error between the exact solution and the appropriate solution we measure the distance between two particles, one transported by the exact flow and the other by the approximate flow. In the process L°° estimates of the finite element solution of the Dirichlet problem for -A will have to be established following the arguments of Nitsche [12]. For the sake of clarity and also because the method of characteristics in the finite element context can be useful for other hyperbolic systems, we begin with a presentation of the method for the transport equation. Then, in Section 2 the method and the error analysis is explained for the two-dimensional Euler equation. Finally, the numerical implementation and some numerical tests are presented in Section 3. 1. Finite Elements and Characteristics for the Transport Equation. 1. a. Statement of the Problem. Let fi be a bounded open set of R", let Q = Q X ]0, T[, and u a divergence free (V-m = 0) vector of (H\Q) n L°°(Q))n. Let/and p0 be two functions of LX(Q) and consider the problem , v ^+ mVp = / in fi x 10, 7T = Q, (1.1) 9r H J J ' L *' p(x, t) = p0(x, t) for all (x, r} e 5 = fi X (0} u 2", where T is the boundary of fi, v is outward normal and 2~ = {{x, t): u(x, t) ■ v < 0, x e T}. 2" represents the part of the boundary of Q where the flow enters into fi; thus License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use THE VORTEX METHOD WITH FINITE ELEMENTS 121 the boundary conditions for p are given at initial time and when the velocity u enters into fi. Several physical phenomena are governed by this equation, known as the transport equation. On rectangular domains fi (1.1) is easily discretized by any upwind finite difference scheme, but if fi is complicated there is no simple nondissipative finite element scheme, and the method of characteristics is usually considered as an expensive numerical method. Let us show that for first order accurate discretizations this is not so costly. The method of characteristics is based upon the following observation: given {x, /)EÖ define [Xx''(t), t} by dX \u(X(r),r) if*(r)efi, (1.2) ^=l0 ifX^fi, Vre]0,r[, X(t) = x. If u is uniformly Lipschitz continuous with respect to x, (1.2) has a unique solution on ]0, r[; then we define (1.3) p(x, t) = Po(X*'(0), 0) + f7(X*-'(t),t) dr and claim that p is a solution of (1.1). If the data u, p0 and/ are not smooth the proof is difficult [15], but if p0,/and u are in C'(£?)> men p is differentiable with respect to {x, r} and p(Xx''(r), t) p(x, 0 "§(t 0 + VxP(X*''(t) x) (1.4) + o(t t) + o(X x) dp and also -(^ + «vj^t o + (X*>'(r), t) p(x, i) = Cf(Xx'(a), a) da. Jt Therefore we have the following result: Proposition 1. When u is uniformly Lipschitz continuous in x and divergence free and p0, m and f are in LX(Q), then the solution o/ (1.1) is given by (1.2), (1.3). Remark. Equation (1.2) includes the possibility that a characteristic leaves fi before t = 0 with the convention that/|r = 0. l.b. Discretization. To discretize (1.1) we choose a triangulation <5h of fi made of nonoverlapping triangles if n = 2 or tetrahedra if n = 3 with the usual properties [4]: \ = kir-, (1 5) ri n ry = 0 or one vertex, or one side (or face), (or one edge), N U t, = Q,h c fi, distance (fiA, fi) = 0(h2), h = size of largest side of %. i License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 122 CLAUDE BARDOS, MICHEL BERCOVIER AND OLIVIER PIRONNEAU If Ar denotes a time step then Q is approximated by (1.6) Q„ = U Pv; Py = T, x ]jAt, (j + l)Ar[, j 1.M £(7/Ar). Instead of u we shall approximate the stream function * of u, i.e. the function such that _/9*3 9*2 9*, 9*3 9*2 _ 9*, \ y 9x2 9x3 ' 9x3 9x, ' 9x, 9x2 / (L7) =[— — ) if« = 2 \ 9x2' 9X| / Let *A be an approximation of * in the space Hh: (1.8) Hh = [rjpA: T) iix G Öa> dt 10 otherwise, for all t E ]tJ, tJ + Ar] with X(tJ) = £,y. Set (1.14) g*M/+« = x^ tJ+i = tJ + Ar, m,j 1 (1.15) ph(e(i)J+\ tJ+l) = Ph(tiJ, ñ + 2 /*(*£ tfX«+i 'i0> i where k(i) is such that (1.16) Ê*W+" E Pk(iy+l. 2. For all /' < N such that there are no k(i) = /' with Procedure P(è''°, tJ+\ tJ+x tj-\-uh, -\) compute all the nodes {**,/*}? of the characteristic solution of (1.13) for all r E ]0, tJ+l] with X(tJ+]) = |''°, and set (1.17) trj+l-ln, "h. (1.18) ph(^+\ r>+1) = Poh(X(0)) + ZMxï. <Ä)('« '«+i)i 3. Replace j by/ + 1 and stop if/ > M, else keep only one |* per t, and go back to 1. Remark. The nodes {xj{, t'm) of the characteristics are the vertices of these broken lines. License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use 124 CLAUDE BARDOS, MICHEL BERCOVIER AND OLIVIER PIRONNEAU Procedure F(£, t, At, vh, e). Computes with the method of Algorithm 1 the intersections [xm, tm) with the boundaries dP0 of the elements of the triangulation of Qh, of the characteristic X(t) solution of dX Í on It, t + AtI ife=+l, on [t — At, t[ ife = -l, — = vh(X,t),oxQÚX <2fi„ dt (1.19) Alternatively we may choose to use with X(t) = f Algorithm 3: Computes the solution ph(x, t) of (1.10), (1.11). 0. Choose 9 = T/r, r E N. 1. For /' = 1 to r compute with Algorithm 2 the solution of (1.10), (1.11), over r E [(/' 1)9, i9[ with ph(-, (i 1)0) for initial condition. Comments. For regular fields uh most elements are expected to contain one end point of the characteristics that are computed forward in time. For those elements which have no such characteristic at a given time we pick any point, here |'°, and compute the characteristic that ends at £'°, backward in time; last, note that it is not absolutely necessary to go back till r = 0; the user may pick a time 0 and ph(t) can be computed from ph(t — 6) according to (1.18); this gives Algorithm 3. These schemes could be made conservative by putting appropriate weights in (1.14), (1.15) according to the number of characteristics that end in the same element and start from the same element, but the error analysis below would no longer be valid. The number of operations, NX, of Algorithm 3 is of the order of I D2T I u.At \ uh9\ (1.20) ^<(^-max(^-,l) + ^)c, where C, is some constant of the order of 10, v is an upper bound on the number of backward characteristics computed at each time step, D is the diameter of fiA. We have found that the best choice for Ar is the one that makes the characteristics cross one element, i.e.,